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ABSTRACT 

The ability of primordial gas to cool in proto-galactic haloes exposed to Lyman- 
Werner (LW) radiation is critically dependent on the self-shielding of H2. We perform 
radiative transfer calculations of LW line photons, post-processing outputs from three- 
dimensional adaptive mesh refinement (AMR) simulations of haloes with Tvh- > 10^ 
K at z ~ 10. We calculate the optically thick photodissociation rate numerically, 
including the effects of density, temperature, and velocity gradients in the gas, as well 
as line overlap and shielding of H2 by HI, over a large number of sight-lines. In low- 
density regions (n < 10'' cm~'^) the dissociation rates exceed those obtained using most 
previous approximations by more than an order of magnitude; the correction is smaller 
at higher densities. We trace the or igin of the deviations prim arily to inaccuracies of (i) 
the most common fitting formula ()Draine fc Bertoldil[l996h for the suppression of the 
dissociation rate and (ii) estimates for the effective shielding column density from local 
properties of the gas. The combined effects of gas temperature and velocity gradients 
are comparatively less important, typically altering the spherically averaged rate only 
by a factor of < two. We present a simple modification to the DB96 fitting formula for 
the optically thick rate which improves agreement with our numerical results to within 
~ 15 per cent, and can be adopted in future simulations. We find that estimates for 
the effective shielding column can be improved by using the local Sobolev length. Our 
correction to the H2 self-shielding reduces the critical LW flux to suppress H2 cooling 
in Tvir > 10^ K haloes by an order of magnitude; this increases the number of such 
haloes in which supermassive (M ~ 10^ M0) black holes may have formed. 
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1 INTRODUCTION 

It has long been known that molecular hydrogen, the most 
efficient coolant in metal-free gas at temperatures below 
lO'^K, played a key role in formation of first-ge neration, 
"Population III," stars fsee I Abel fc Haimanll2000l for a re- 
view). As soon as these first stars began to shine, however, 
they also began to destroy H2 via dissociating (LW) pho- 
tons in the range 11-13.6 eV, to which the universe is largely 
transparent even at early times, z ^ 20 — 30. The nature and 
extent of this photodissociation feedback has important con- 
sequences for subsequent star-formation, reionization, and 
the formation of massive black holes at early times. 

In regions where large H2 column densities build up 
(A'^H2 > 10^* cm~^), photodissociation is suppressed as 
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the LW bands become optically thick; the cooling proper- 
ties of UV-irradiated primordial gas thus depend largely 
on its ability to "self-shield." Unfortunately, the prob- 
lem of modeling self-shielding exactly in existing stud- 
ies remains intractable; in particular, approximate treat- 
ments are necessitated by two main challenges, which 
are the focus of this paper. First, the computational ex- 
pense for three-dimensional simulations of finding the ex- 
act self-shielding column density in a large number of 
directions is prohibitive. As a result, studies have typ- 
ically either adopted the op t ically - thin d issociation rate 
throughout llMachacek et al.1 200 ll. l2003l: iMesinger et all 
'20061. 12OO9I ; I Wise fc AbellbOO?! . l2008al lbl; iGreif et al.ll2010l ). 
relie d upon estimating A^H9 from local prope r ties of the gas 
Bro mm fc Loebl |2003|; I Johnson et al.1 Booj; iGnedin et al.1 



120091 ; IShang et all I2OI0I ; I Johnson et al.1 l201ll ). or sacri- 
ficed angular resolution, finding the exact column density 
in a small number of directions to estimate the dissoci- 
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ation rate (lYoshida et all |2003| . l2007l : iGlover fc Mac Lowl 
l2007al lb[). Alternatively, some have employed a local method, 
which allows for contributions to shielding only from 
gas within a single smooth-particle- hy drodynamics (SPH) 
smoothing length (|Glover et alj |2006| ) or within a width 
defined by the size of the underlying simulatio n grid - a 
method also investigated by iGlover fc Mac Lowl (l2007al lbh. 
Recently, an algorithm for finding the projected column 
density distribution as seen by each SPH particle, us- 
ing a Healpix tessella tion with 48 equal-ar ea pixels, has 
been implemented by IClover fc ClarkI l|201lll . In one-zone 
models, A^H2 must be specified fr om only "local" prop- 
erties of the gas, by definition 



2 NUMERICAL METHOD 



2.1 Simulations 
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Even once an estimate of the self-shielding column is 
obtained however, finding the exact photodissociation rate 
represents a large computational expense, requiring high nu- 
merical resolution in order to explicitly account for pro- 
cessing of the incident LW radiation as a function of fre- 
quency. Furthermore, even when only the (ortho and para) 
ground states of the molecule are populated, there are al- 
ready a total of 76 LW transitions that contribute to the 
total optical depth in the relevant frequency range (pho- 
ton energies below 13.6eV). Most often, studies circumvent 
this difficulty by a d optin g analytic expressions provided by 
iDraine fc Bertoldil (Il996l . hereafter DB96) to model self- 
shie lding, with a few exceptions among semi-analytic mod- 
els jHaiman et al.ll2000l : ICiardi et al]|2000l : IGlover fc Brand 
2001, 200^ and one-dimensional simulations (^Ricotti et al 



2001 : Hosokawa fc Inutsukal2005l . lioOS ') , which include more 



detailed calculations. 

The two primary goals of this study are (i) to quan- 
tify the accuracy of previous models for self-shielding, and 
(ii) to provide an improved analytic fit for the suppres- 
sion of the photodissociation rate by shielding, which can 
be used in future simulations. To accomplish this, we post- 
pr ocess the out p uts fr om a suite of simulations performed 
bv IShang et all (|2010l . hereafter SBHIO), who studied the 
effects of UV-irradiation on protogalactic haloes with virial 
temperatures Tvir > 10*K. Within these haloes, we calculate 
the exact optically-thick dissociation rate at a large num- 
ber of points, in a large number of directions, with a detailed 
treatment of radiative transfer, explicitly including the ef- 
fects of density, temperature, and velocity gradients in the 
gas, which are most often neglected in existing models. 

The rest of this paper is organized as follows. In § [2l we 
describe the SBHIO simulations and our methods to com- 
pute photodissociation rates with shielding. In this section, 
we also recapitulate several of the most common existing 
methods for modeling self-shielding. In § |3l we present our 
main results - the numerically computed suppression of the 
photodissociation rates. We compare these results to those 
obtained with each of the previously used methods, and we 
also elucidate the effects of gas temperature and velocity 
gradients, as well as the accuracy of the analytic expres- 
sions from DB96. In § IH we discuss the implications of our 
results and the associated uncertainties, and we offer our 
conclusions in § [5] 



We utilize outputs from a suite of simulations performed 
by SBHIO with the Eulerian adaptive mesh refinement -I- 
N-body code ENZO (jBrvanI [l999l : iNorman fc BrvanI Il999l : 
lO'Shea et al.l 12003 ). For the sake of brevity, we hmit the 
discussion here to the most pertinent features of their nu- 
merical method, and refer the reader to the original study 
for further details. 

The simulations were performed within a comoving box 
1 Mpc on a side and assuming a ACDM cosmologi- 
cal model with standard concordance parameters: S7dm ~ 
0.233, Qi, = 0.0462, = 0.721, as = 0.817, n, = 

0.96, and h — 0.701. A preliminary run was initialized with 
a root grid of 128"' and no nested grids. This was performed 
(with radiative cooling turned ofi^) in order to identify haloes 
with virial masses of a few x 1O^M0 at 2 ~ 10. Three of 
these halo were then re-simulated at high resolution with 
new initial conditions - three nested grids with an efi'ective 
innermost resolution of 1024"^ - and with levels of refine- 
ment added adaptively. Refinement was increased when the 
baryon or dark matter mass exceeded thresholds of 68 and 
683 Mq respectively, and in order to maintain sufficient res- 
olution of the local Jeans length - at least four grid cells - to 
prevent artificial fragmentation. Throughout, the dark mat- 
ter (DM) gravity was smoothed on a scale of 0.954 (co- 
moving) parsec. Each simulation was stopped after reaching 
a refinement level of 18 - a resolution of 0.0298 parsec 
(comoving) or 800 AU (physical). 

The non-equilibrium chemistry for a gas of primordial 
composition was followed with a chemical network com- 
prising 28 gas-phase reactions, including H2 photodissoci- 
ation. Radiative coolin g by H2 was modeled with the func- 
tion provided bv iGalli fc PaUal |l998 ). Each halo was run 
with varied LW backgrounds in a range of intensities J21 — 
1—10^, where the standard normalization is implied here and 
throughout, Jlw ~ J21 x lO"'^^ erg s~^ cm~^ sr~^ Hz~^, 
and Jlw is the fiux intensity at the average LW frequency 
{hv = 12.4 eV). The photo-dissociating source was modeled 
as a blackbody with a temperature of 10* or 10^ K (referred 
to as T4 and T5 respectively). Since in the T4 case, the 
photodissociation of H~ , rather than direct H2 dissociation, 
controls the H2 chemistry, we restrict our analysis here to 
the T5 runs. We defer further details of the adopted self- 
shielding model to § 12.31 



2.2 Numerical self— shielding calculations in 3D 

From the simulation outputs, we selected five snapshots 
spanning a range in redshift z ~ 8 — 12, in which the gas is 
optically thick to LW photons. Three of these - outputs 1, 
2, and 5 - are designated "cold," as they were subjected to 
a modest dissociating fiux, and thus have experienced sig- 
nificant H2-cooling. In outputs 3 and 4 the gas remains hot, 
with temperatures > 7000 K, having been exposed to a very 
strong LW flux that kept the gas H2-poor. The physical 
properties of the five halo snapshots are summarized in Ta- 
ble [T] and the radial profiles of the density, temperature, H2 
fraction, and the electron fraction are shown for each snap- 
shot in Figure[T] We selected ~ 100 points from the cold halo 
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Figure 1. Spherically averaged profiles of temperature, H2 frac- 
tion, electron fraction, and particle density in the selected simu- 
lation outputs. The radius, R, is measured from the densest point 
in the halo. 



Table 1. Physical properties of haloes selected from the SBHIO 
suite of simulations. Throughout, the outputs will be referenced 
by number; 1-4 correspond to "Halo A" in SBHIO, and #5 to 
"Halo C." The redshift {z) of each output is given, along with 
the virial mass at the collapse redshift (mvir coi; a-s in SBHIO, the 
virial mass is defined as the total mass in baryons and dark matter 
within a spherically averaged overdensity of 200 with respect to 
the critical density). The temperature and total particle densities. 
To and no, respectively, are specified at the densest point in the 
halo. The intensity of incident Lyman-Werner radiation, Jlw, in 
each run is parametrized in the standard manner: Jlw = 1^21 x 



10-21 erg 
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cm 2 Hz 


1 sr- 


1 






output 


J21 
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nivir 


col(M0) 


ro(K) 


no (cm 


1 


10^ 


12.274 


2.45 


xlO' 


9.2 X 10^ 


3.1 X 10** 
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10^ 


12.285 


2.45 


xlO'^ 


4.0 X 10^ 


3 X 10"' 


3 


10^ 


9.93 


5.49 


xlO'^ 


6.3 X 10^ 


1.4 X 10^ 


4 


10^ 


9.94 


5.49 


xlO'^ 


7.8 X 10^ 


1.9 X lO"* 


5 


102 


8.3 


7.89 


xlO'^ 


9.3 X 10^ 


1.6 X 10** 



outputs at which to calculate the three-dimensional H2- 
dissociation rate. Fewer points 30) from the "hot halo" 
outputs were analyzed, as these are considerably more ho- 
mogeneous with respect to temperature, density, and chem- 
ical composition (and are nearly optically thin even in the 
most dense regions). The selected points span a range in 
radii from 0.1 — 10 pc (physical) where the radius is defined, 
here and throughout, as the distance from the densest point 
in the halo. The number density and temperature at the se- 
lected points vary between 10 — 10 cm" -3 and 300 - 10*K, 
respectively. 



cited from the electronic ground state, X^E^, to the B^Ej 
or C^Ilu stat^B) the Lyman and Werner bands respectively. 
Subsequent decays lead to the vibrational continuum of the 
ground state, rather than to a bound state ~ 15 per cent of 
the time, thus dissociating the nuclei. 

The "pumping rate" from a given rovibrational state 
(i;, J) to the excited electronic state with (v' , J') is: 



Jv 



(1) 



where is the frequency dependent cross-section and ftp 
is Planck's constant. The frequency threshold, Uth, corre- 
sponds to the lowest energy photons capable of efficiently 
dissociating H2, with hf « 11.1 eV. We do not include flux 
at hu ^ 13.6 eV, as photons with energies above the Ly- 
man limit are assumed to have already been absorbed by 
the neutral HI in the intergalactic medium (IGM) outside 
the halo (and were not included in SBHIO). 

The dissociation rate from the initial {v, J) is then ob- 
tained from the product of the pumping rate and the fraction 
of decays leading to the vibrational continuum from («', J'), 
with a sum taken over all possible upper states: 



^ ^ Cv ,J ,v' ,J' fdi 



iss,u' , J' • 



(2) 



The dissociation p r obabi lities here, /diss,i;', j' , are obtained 
from lAbgrall et al.l l|2000l ). The sum over rates from all lower 
levels, weighted by the fraction of molecules initially in each, 
fv,j, then gives the total rate: 



fed 



— ^ ^ ^diss.v ,jfv 



(3) 



where the fv.j are given according to a Boltzmann distribu- 
tion, unless otherwise specified. 



2.2.2 Radiative transfer in the haloes 

In the three-dimensional calculations of fcdiss, the first step 
is to generate a set of rayfl emanating from each point where 
the dissociation rate is to be found. These sample evenly in 
the azimuthal angle and in the cosine of the polar angle, 
tiling a sphere of radius ~ 100 pc. Along each ray the prop- 
erties of the gas are sampled at intervals determined by the 
size of the underlying grids; the distance between sample 
points is typically 0.01 pc in regions with the highest levels 
of refinement, and increases to ~ 5 pc toward the outskirts 
(areas with lowest resolution) of the halo. 

In our fiducial calculations, the spectrum of the inci- 
dent radiation is initialized to be flat in the range 11.1-13.6 
eV, with the implicit assumption that processing of the LW 
background in the IGM is negligible. Note that, in general, 
the cosmological background will be modifle d, though pri- 



maril y by HI absorption, and no t by H2 itself 



i2000l , hereafter HAROO, see also lRicotti et all 



Haiman et al.l 



20011 ) and the 



impact of the resulting "sawtooth" modulation is considered 



2.2.1 Rate of photodissociation 

H2 is photodissociated primarily via th e two-step S o lomon 
process (Solomon 1965 ; see also iField et ah I 1 19661 : 
IStecher fc WiUiamd 1 19671 ). in which molecules are ex- 



1 C^Tlu is split into HJ and H^ states owing to A-doubling, the 
two-fold degeneracy of each rotati onal level (J). 

We use the analysis toolkit YT l lTurk et al.ll201(]|) to interface 
with the raw simulation data; see § 13.51 for the required angular 
resolution. 
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in ij |4.2l Tracing a ray from tiie outside in, the gas is treated 
as a series of thin slabs, each with uniform density, tempera- 
ture, bulk velocity, and chemical composition defined at the 
sample point; the column density is specified by nn^ x As, 
where As is the width of the slab. The frequency-dependent 
optical depth, r,^, of the slab is then obtained by summing 
over contributions from all included LW transitions, each 
of which is modeled by a Voigt profile. The rest-frame fre- 
quencies are Doppler shifted according to the slab's line- 
of-sight velocity relative to the point where the dissociation 
rate is to be calculated. The numerical wavelength resolution 
(AA = 2 X 10~*A at the lowest temperatures) is set adap- 
tively in order to always resolve the thermal line width and 
is sufficient to explicitly account for overlap of the Lorentz 
wings. The nec essary molecul a r data for these calculations 
are provided bv lAbgrall l|l993al lb[). 

We include transitions from the 29 bound rotational 
levels within the ground electronic and vibrational {v — 0) 
states to excited states with v' ^ 37, J ^5 10; in total, this 
amounts to 1492 possible transitions out of i; = 0. We do 
not include absorption from higher vibrational levels, which 
are populated only at particle densities much larger th an we 
consider her e, n > 10 cm" jLe Bourlot et al.lll999l . here- 
after LPF99. [Flower fc Harrisll2007l 'l. In our fiducial calcula- 
tions, the rotational levels of u = are populated according 
to a Boltzmann distribution defined by the temperature of 
the slab, irrespective of the local density. Strictly speaking, 
populations in rotational states within w = do not ther- 
malize until (temperature dependent) critical densities are 
reached, usually taken to be nc rit — 10* cm~"^ (but see Ta- 
ble 1 in iFlower fc Harris! [2007I ). at which depopulation of 
excited states is dominated by coUisional de-excitation|f| To 
further address this issue, we perform an additional set of 
calculations in which all molecules are assumed to be in the 
ground states of para (ortho) hydrogen, v = 0, J = (1). In 
this case, all 76 possible transitions are included; the results 
are discussed in S 14.41 

Finally, having found along a sightline, the photodis- 
sociation rate is calculated at the point of interest. This pro- 
cedure is repeated for each of the sight lines and the final 
rate is obtained from a simple average over all directions. 
Note, however, that the spherical average is only meaning- 
ful in the case of an isotropic UV background. This would 
be relevant when the Olber's integral for the local flux is 
dominated by a large number of distant sources. For haloes 
with unusually bright and close neighbors, the flux can be 
dominated by a single (or a few) of the nearby sources. In 
ij |4.3l we discuss a radiation field with a preferential direction 
that would be relevant in this case, and interpret our results 
in this context. 

One remaining caveat is the possibility of H2 shielding 
by HI. In general, Lyman-series absorption within the haloes 
can suppress the H2 dissociation rate by a large factor, 
though this requires high optical dep th in the wings of the HI 
lines, and thus A'^hi > 10^^ cm~^ (^Wolcott-Green fc HaimanI 
[20l3). Sufficiently large neutral column densities are indeed 
present in the the outputs we analyze, though only at small 



^ Note that the critical density is different for each species that 
perturbs the molecule; here and throughout, we refer to n^rit for 
collisions with atomic hydrogen only. 



radii (< a few pc), and the resulting modification of fcdiss in 
these regions is discussed in § 14.21 

2.3 Approximate treatments of self-shielding 

Approximations for self-shielding in existing simulations 
have been necessitated by the two challenges mentioned 
above: first, the computational expense for simulations of 
finding the H2 column density in a large number of direc- 
tions is prohibitive. An estimate for this inherently non-local 
quantity therefore typically must be obtained from purely 
local information. In § 12.3.11 we give a detailed account of 
several ways in which this is commonly achieved. Second, 
calculating the exact suppression of the optically-thin rate, 
with full radiative transfer in each LW line, is expensive, 
requiring high numerical wavelength resolution, as well as 
the inclusion of non-local effects from the temperature and 
velocity structure in the gas. Therefore, studies often rely on 
an analytic expression for the optically-thick rate provided 
by DB96. We briefly describe this method in ij |2.3.2l includ- 
ing the assumptions and limitations in applicability of the 
analytic fitting formula. 

2.3.1 The self-shielding column density 

In order to estimate the column density, several common 
methods make use of local properties of the gas to define a 
characteristic length scale, Lchar; the column density is then 
obtained from: 

A*'h2 = f^H2-^'char, (4) 

with the implicit assumption that the H2 number density, 
riHa, is constant at the local value over the length Lchar. 
Several methods have been used to define a characteristic 
length scale: 

The Jeans Length 

A common approach in both simulations and one-zone 
models is to assume the total mass in the optically-thick 
region is of order the Jeans mass, with ichar then defined 
using the local Jeans length. In regions where particle den- 
sities exceed ntot > 10^ cm"'', SBHIO have shown that this 
provides a very accurate estimate of Nu2', however, in lower 
density regions, it typically overestimates the integrated 
by up to an order of magnitude (see their Figure 9). 
This method also does not account for temperature or 
velocity gradients in the gas. 

The Sobolev Length 

Large velocity gradients will cause the LW resonances, as 
seen from a given point in the cloud, to appear shifted from 
their rest-frame wavelength, making the systematic deple- 
tion of dissociating photons less efficient than if the gas were 
static. If the fluid motions are largely disordered and/or su- 
personic, the exact column density along a line of sight can 
then exceed the effective self-shielding column density, A^off, 
by a large factor. 

In general, A'efr depends on the detailed velocity struc- 
ture in the gas; however, it may be estimated by taking 
a characterist ic length equal to the Sobolev length, Lsob 
l|Sobolevlll95"7l ). This defines a distance over which the mean 
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(macroscopic) fluid velocity changes by a factor of the (mi- 
croscopic) thermal velocity of the molecules, vth, 

r — 

|dii/ds| 

with the assumption that dt;/ds is constant. Then for a given 
sightline, Lsob is the distance at which the absorption pro- 
files will be shifted by one local line width. In adopting 
the Sobolev length for L^har, it is implicitly assumed that 
Doppler shifts in the LW resonances of molecules closer than 
Lsoh can be ignored, and that all molecules beyond this dis- 
tance contribute negligibly to shielding. 

Strictly speaking, this method should only provide an 
accurate estimate of A^eff if the velocity gradient is both large 
and monotonic and the intrinsic line widths are negligible 
in comparison to the Doppler cores. If instead the Lorentz 
wings (and line overlap) become important, these will render 
the Doppler shifts irrelevant and this method will underes- 
timate self-shielding as a result. None the less, it may prove 
useful, particularly if motions of the fluid are supersonically 
turbulent. In fact, a similar method has been fruitfully em- 
ployed in the analogous case of t he escape fraction of p hotons 
from dense, metal-free gas (e.g. lYoshida et .aLll2006l ). 

In order to generalize to a non-spherically symmetric 
geometry, one option is to define a single Lsob from the 
mean Sobolev length over all directions, which we will refer 
to as Lgob- Alternatively, in direct analogy with the three- 
dimensional self-shielding calculation, we could find Lsob 
and the corresponding dissociation rate for each sightline 
and take the mean rate over all directions (hereafter Lsob)- 
Clearly, this approach should provide a reasonably accurate 
estimate of the dissociation rate if the corresponding column 
density is a good approximation for A^off along every direc- 
tion. Finally, we could take the minimum Sobolev length 
(Lgjj'b), or an average column density (isob)) over all direc- 
tions. 

A "Sobolev- Like" Length 

A method akin to the Sobolev length has been recently em- 
ployed and shown to provide a fairly accurate estim ate of the 
integrated column density by iGnedin et al.l (|2009l '). though 
in a different context than that of the present work (they 
model self-shielding in individual star-forming regions of a 
Milky Way progenitor). In this case, a characteristic length 
is obtained from 



LL 



|Vp|' 



(6) 



thereby defining a distance over which the gas density, p, 
should be significantly diminished, and assuming that the 
optical depth beyond isob is negligible. 

We expect this method to prove most useful when 
the length scale of significant decrease in the gas density 
is shorter than those over which large variations in the 
H2 fraction or the gas velocity occur. Note that the aim 
is to define a distance beyond which the gas contributes 
negligibly to shielding, so extrapolation along lines of sight 
for which dp /As > is not meaningful (this happens for 
off-centre points, in directions toward the halo centre). 
Thus, we define I/gobi -^'Sobi S'Hd ^s^b in ^ similar manner 
as described above for the traditional Sobolev-like length, 
but include in the averaging only those sight-lines for which 
dp/ds < . 



The "Six-Ray Approxtmation" 

Finally we consider one non-local method, variations 
of whi ch have b e en im plemented i n simulations by, 
e.g.. lYoshida et all (|2003l . HoOT,) and iGlover fc Mac Lowl 
(|2007al lbl). In this case, the exact is obtained by in- 

tegrating the H2 profile along six lines of sight parallel to 
the Cartesian axes. The value of the shield factor for each is 
obtained from one of the DB96 fits (equations [8] or |9] below) 
and the final rate is found by averaging over six directions|j 
In spite of the low angular resolution, this approach is ex- 
pected to be reasonably accurate unless the gas is (super- 
sonically) turbulent, in which case neglecting Doppler shifts 
of the LW lines will likely cause the integrated A^h2 to sub- 
stantially exceed the effective self-shielding column density. 



2.3.2 Analytic approximations for the photodissociation 
rate 

Self-shielding by H2 has received attention over the years 
for its importance in the context of interstellar clouds (e.g . 
HoUenbach et alj Il97ll: IShullI Il978l: iFederman et al l 1 1979^ : 



de Jong et al.l Il980l: lAbgrall et al.Ml992l : iHeck et al.l Il992l : 

Le Bourlot et al]|l993l : DB96). A number of analytic models 
for the attenuation of the incident flux have been put for- 
ward, as explicit calculations of the radiation fleld as a func- 
tion of both cloud depth and frequency have not previously 
been feasible. In the last decade, however, studies have most 
often employed the expressions provided by DB96. These au- 
thors model a semi-infinite, static slab of gas irradiated on 
one surface, and parametrize the dissociation rate with a 
"shield factor" as: 

fcdis.(iVH2,T) = /.h(iVH2,r) X fcdi..(iVH2 = o,r), (7) 

where A;diss(A'H2 = 0) is the optically-thin rate. They show 
that, at low temperatures (T ~ a few xlO^ K), suppression 
of the optically-thin rate can be well approximated by a sim- 
ple power-law that depends only on the H2 column density 
as: 



/sh,36 (A*'h2 



-3/4 



(8) 



These authors also provide a slightly more complicated func- 
tional form, which attempts to incorporate a temperature 
dependence due to thermal broadening of the lines, and 
which fits their results more accurately: 



/sh.37 {N^„T) = 



0.965 



0.035 



X exp [-8.5 X 10"'' (1 + xf-^] . (9) 

Here x = NH2/5 x 10^* cm"^, 65 = 6/IO'' cm s"\ and b is 
the Doppler broadening parameter (equations [S] and d] are 
given in DB96 as their equations 36 and 37, respectively, as 
indicated by the subscripts). Because]^ is the more accurate 
of the two, this will be the focus of our discussion henceforth. 

These expressions have been ubiquitously used to model 
self-shielding in simulations. While it is often noted that 



* Note that lYoshida et al. I 1I2OO3I ') use instead the shield factor 
calculated for the line of sight with minimum H2 column density. 
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Figure 2. The numerical results for the shield factor in three-dimensions, /sh,3Di are shown at the local temperature (left-hand panel) 
and particle density (right-hand panel). Points are included from all five halo outputs; those from cold and hot haloes are denoted (here 
and in all subsequent figures) by open triangles and squares, respectively. 



they are only strictly valid in the static limit, it should be 
emphasized that DB96 make several assumptions in their 
modeling which limits the apphcability of the expressions 
they derived. Most importantly, they consider a cold gas, 
with T < a few x 10^ K. Furthermore, they assume either 
(i) a steady-state rovibrational distribution that includes 
absorptions from u > or (ii) an isothermal gas with H2 
level populations given by a Boltzmann distribution. Also 
included are shielding by dust and the formation of H2 on 
grains. None the less, their fitting formulae are often imple- 
mented in contexts quite different from the original study, in 
which these assumptions are not well motivated. In hght of 
this, §[3]addresses the accuracy of equation ([9]) in the present 
context, i.e. for metal-free gas with a wide range of densities 
10 < ntot < W cm~^ and temperatures 300 < T < 10*K, 
and in which dust, as well as UV-pumping from excited 
vibrational states is likely negligible (see 14.41) . 



3 RESULTS 

The numerical results for the photodissociation rate in 
three-dimensions, parametrized (as above) by a shield fac- 
tor, /ah, 3D, are shown in Figure [2] as a function of the local 
temperature and particle density of each point. The trend 
toward decreased optical depth at high temperatures and 
low densities is due to the characteristic structure in these 
haloes, in which regions farther from the dense core are both 
hotter and more rarefied, and typically see lower shielding 
columns (see Figure [T]). The exception is the hot haloes, in 
which the gas remains optically thin even at small radii, 
where n > 10"^~* cm~^, due to the low molecular fraction. 
As a result, these hot haloes (outputs 3 and 4) will not be 
further addressed in great detail, as they add little to the 
discussion of self-shielding. 

Figure [3] shows these results in comparison with several 
approximate methods. The four panels show /^h given by 
equation ((9|, and in each, the column density is specified by 
one of the methods discussed in ii l2.3.1l Because studies typ- 
ically must employ estimates for A'h2 as well as an analytic 
fit for the shield factor, the overall discrepancies seen here 
are representative of what would be realized in a simulation. 

A clear trend is immediately apparent in Figure [3l 



which appears independent of how the column density is es- 
timated: all approximate methods significantly overestimate 
shielding compared to our numerical results in the range 
10"'^ < /sh,3D < 0.3, while each is considerably more ac- 
curate at smaller values of /ah, 3D . In principle, this discrep- 
ancy may arise from physical effects, i.e. Doppler shifts in 
the LW absorption lines, causing self-shielding to be weaker 
than in a static gas, or variations in the gas temperature, 
resulting in depopulation of states that contribute most to 
dissociation. Alternatively, it may be due to inaccuracy of 
the analytic fit for the /ah itself. In what follows, we examine 
each of these possibilities in turn, to determine how much 
each effect contributes to the observed trend. 



3.1 Analytic approximations for /sh 

The left panel of Figure [4] shows a comparison of the shield 
factor obtained from the DB96 expression (equation [S]), 
against the results from our numerical calculations at T = 
500, 1000, and 5000K. In the latter, the gas is modeled as a 
static and isothermal slab, in order to isolate the accuracy 
of equation ((9| from the effects of temperature and velocity 
gradients. It is apparent from this figure that the trend of 
discrepancies seen in all four panels of Figure [3] is caused 
in large part by the inaccuracy of the fitting formula itself 
at temperatures greater than a few hundred Kelvin. (Note 
that, as shown in Figure [21 all points in the simulation with 
/ah, 3D > 10""^, which are discrepant, are at T > 500 K.) 

We have found that the large discrepancies - equa- 
tion underestimates the numerical results by up to an 
order of magnitude - are due to the temperature dependence 
of /ah for thermalized H2 populations. To illustrate this, con- 
sider a gas in local thermodynamic equilibrium (LTE) at a 
(uniform) temperature of a few hundred Kelvin, so that only 
the lowest rotational states within the vibrational ground 
state will be significantly populated (for reference, the en- 
ergies of the « = 0, J = 1,2 states are « 170,509 K). As 
the temperature is increased, the populations will be diluted 
over a greater number of rotational levels (e.g. at T ~ sev- 
eral thousand Kelvin, non-negligible populations build up 
in J < 15). The upshot of thus spreading absorbers over a 
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sh,3D ^ sh,3D 



Figure 3. Comparison of our numerical results for the average shield factor in three dimensions, /sh,3Di with those given by the analytic 
fitting formula (/gjj 37 in equation l[9}, taken from DB96). Clockwrise from the upper left, the column density is specified by the best— fitting 
Sobolev and Sobolev-like ( "Sobolev-p" ) methods, (ig^j^ and L^^^ respectively; see 12.3. Il l, and the Jeans length. The lower left panel 
shows results of the only previous non-local approach, the "six-ray method," in which equation JgJ is used to calculate the shield factor 
along six sight lines, and an average over these gives /gh 37(six-ray). Non-local effects of velocity and temperature gradients are included 
in the calculation of /sh.3D 1 described in § 12.21 



greater number of states is that shielding becomes weakei[f|. 
This effect is not modeled by the DB96 expression (equa- 
tion [§]), which accounts only for decreased shielding due to 
thermal broadening of the lines. As a result, we find that it is 
accurate only for a rotationally-cold gas, i.e. when molecules 
occupy only the first few J states (we show this explicitly in 
Figure [TT|) . 

The inaccuracy of this ana lytic fit at high temper atures 
has been noted previously by lAhn fc Shapirg 112003), who 
proposed that it may be remedied by artificially increasing 
the thermal broadening parameter used in equation ((9|. We 
have found that this method yields very little improvement 
in accuracy; however, a much better fit to our numerical 
results is obtained with only a slight modification of equa- 
tion ^ as follows. We treat a, below, as a free parameter, 
where a = 2 in the original expression: 



/sh (iVH,,r) 



0.965 



0.035 



X exp [-8.5 X lO""* (1 + a;)°-^] . (10) 
We find that a — 1.1 improves the fit drastically in the 

^ The importance of diluti ng the populations over higher energy 
states was pointed out by iFederman et al] l|l979t l , who showed 
that for a thermal distribution, including absorption from J = 2 
changed their results by ~ 40 per cent compared to those in ortho 
and para ground states were populated. 



high temperature regime, with little accuracy sacrificed at 
low temperatures. Physically, this result makes sense, since 
reducing a weakens the temperature effect, which, as we 
argued, is overestimated in the original equation. The modi- 
fied expression agrees with the numerical results to within a 
factor of two at 500 < T < 5000 K, and Nh_2 10^° cm"^, 
as shown by the solid curves in Figure |4] Furthermore, as 
the figure shows, the largest discrepancies occur (a) at low 
temperature, T = 500 K, and larger values of the shield 
factor, fsh > 10~^ or (b) high temperature, T > 5000 K 
and fsh < 10~^. In the context of our study, these combi- 
nations are, in fact, never physically realized, as low tem- 
peratures (T ~ 500 K) are only reached via H2 radiative 
cooling, implying that in cold gas, large column densities 
will have built up and the shield factor will be small; for 
the same same reason, the gas cannot remain hot once it 
has become strongly self-shielding (see the temperature and 
density profiles in Figure [T|. Overall, we find that the agree- 
ment between equation (|10p . with a = 1.1, and the results 
for /sh,3D in the haloes is « 15 per cent. 



3.2 Impact of velocity gradients 

As previously described, the existence of internal velocity 
gradients will cause a given parcel of gas to see more flux 
than it would in a static gas, owing to Doppler shifts of 
the LW resonances in regions of the fluid moving at large 
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Figure 4. The results from equations (|9l left panel) and l|10l with a = 1.1; right panel) are compared against the numerical results for 
/ah at T = 500, 1000, and 5000K. In the latter, the gas is treated as an isothermal and static slab, irradiated uniformly on one side. Tick 
marks along each curve indicate ten— fold decreases (left to right) in the column density, with N-^^ = 10^*^ cm~^ at the lower left; the 
dotted curves are 45-degree lines for reference. 



relative velocities. This can lead the actual A'^h2 to exceed 
the efTective column density by a large factor, and could 
thus also contribute to /sh(A'^approx) overestimating shielding 
compared to the numerical results (see Figure [3]). 

To quantify the magnitude of this effect, we repeated 
our three-dimensional calculations with gas velocities arti- 
ficially set to zero everywhere. The results, shown in the 
left panel of Figure [5] can be divided roughly into three 
regimes: strongest, intermediate, and weakest shielding; the 
transition between the first (latter) two occurring at /sh.so — 
4 X 10~^(2 X 10~^). The dissociation rate is increased in the 
non-static case by an average of 6, 65, and 20 per cent, re- 
spectively, in the three regimes (though in the intermediate 
case, certain points see dissociation rates that are factors of 
up to 4-5 times greater due to relative gas velocities). 

In order to understand why this effect operates differ- 
ently in the three regimes, let us first consider points where 
shielding is weakest, /sh,3D > 2 X 10"^ In Fi gure [6] the 
shield factors are shown for individual sight lines emanat- 
ing from a single point in this regime; these are shown for 
both the static and non-static cases (closed and open tri- 
angles, respectively) at the integrated Nn^ in each direc- 
tion. In the lower panel, the relative line of sight velocities 
along all 16 directions are shown. Three directions are dis- 
tinct from the rest, in that they see the largest velocities as 
well as the greatest integrated column densities - the three 
bottom curves in the lower panel correspond to the three 
right-most points in the upper panel. This highlights an im- 
portant point about the dynamics in these simulations; a co- 
herent flow of the fluid onto the densest region dominates the 
gas velocities, rather than turbulence. Thus, large frequency 
shifts in the H2 lines occur only along sight lines aligned with 
this flow, and while they clearly do change the dissociation 
rate in these directions, variations in these smallest values 
of the shield factor have little effect on the spherically av- 
eraged value. In the strongest shielding regime, the column 
densities are sufficiently high that large damping wings of 
the lines render frequency shifts irrelevant. We have verified 
this by artificially decreasing the column density of each slab 
along these sightlines by a factor of 100 and repeating the 
full calculations with velocities both on and off; in this test, 
/sh,3D becomes several times larger in the non-static case. 



In the intermediate regime, the damping wings of the 
lines are not significant, and frequency shifts have a larger 
impact on the dissociation rate as a result. This is not, how- 
ever, because the total column densities are lower than in 
the previous case. The integrated A''h2 's along sight lines 
emanating from these points are, in fact, comparable to 
those in the strongest shielding regime, but the points them- 
selves are in regions with much higher temperatures (i.e. 
T > lO'^K, compared to a few hundred Kelvin; see Figure 
[2|. This highlights the importance of coupling between the 
effects of temperature and velocity gradients. As previously 
discussed, increased temperature dilutes the II2 populations 
over a greater number of rovibrational levels; the upshot here 
is that diminished Lorentz wings of transitions originating 
from these no longer cancel the effects of Doppler shifts. In 
other words, the column density in a given rovibrational level 
is smaller than in the strongest shielding regime, though the 
total A^H2 for individual sight lines is comparable. 



3.3 Impact of temperature gradients 

While not as often addressed as the effect of frequency 
shifts, variations in the gas temperature along a given line 
of sight may also significantly alter the self-shielding be- 
haviour, particularly for a gas in LTE. There are two dis- 
tinct temperature-induced effects; (i) changes in the level 
populations (specific to the LTE case), and (ii) changes in 
the widths of individual lines. 

Consider the population effect first. For a given parcel 
of gas with temperature To, along a given direction, devia- 
tions in the temperature from To and the resulting change in 
the rotational distribution may either increase or decrease 
the populations of the levels from which most dissociating 
transitions originate. In general, if (a) To is small (~ a few 
X lOOK) and the temperature is increased along a sightline, 
or (b) the difference between the temperature of the shield- 
ing gas and To is very large, the result is to depopulate the 
states from which most UV-pumping occurs, and thus to 
decrease the effective column density. On the other hand, if 
the temperature variations along a sightline are not large, 
the shielding may be increased compared to the isothermal 
case, owing to boosted populations in states from which the 
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Figure 5. Left-hand panel: the fiducial numerical results, /sh,3D, are compared to tliose in which the gas is assumed to be static. Middle 
panel: isolates the efli'ect of setting the temperature constant along each sightline (with the gas again assumed to be static along both 
axes). Right-hand panel: shows the combined eff'ect of the velocity and temperature variations, by comparing the fiducial results with 
those in which the gas is assumed to be both isothermal and static. Note that, here, "isothermal" does not imply that temperature is the 
same for all points, but rather that T = Tq along all sight lines, where Tq is the temperature where the dissociation rate is calculated. 
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Figure 6. Upper panel: comparison of the numerically-calculated 
shield factors along individual sight lines from a single point 
(triangles) with those obtained when the gas is assumed to be 
static (diamonds) or both static and isothermal (squares). The 
full spherically averaged photodissociation rate is increased 23 per 
cent (30 per cent) compared to the static (static and isothermal) 
case. Lower panel: relative line of sight velocities along the same 
directions as shown in the upper panel. The three lower (solid) 
curves correspond to the three sight lines with largest total AfH2 
(closed symbols) in the upper panel. 



strongest transitions (i.e. those for which the product of the 
oscillator factor and dissociation probability is largest) orig- 
inate. 

Interestingly, we find that the line-width effect may 
be equally or more important than the above population- 
induced changes. Consider again our gas parcel with temper- 
ature To, shielded, for simplicity, by a single slab of gas at 
T < Tq. The narrowed thermal widths of absorption lines in 
the slab result in overall much weaker shielding, in compar- 
ison to the isothermal case (though the opacity at the line 
centre is actually increased). Furthermore, because this ef- 
fect may act in the same direction as the population-induced 



changes, shielding can be greatly reduced for the case in 
which T < To- This is illustrated in Figure |6l in which the 
open squares show the results when both the velocity and 
temperature are held constant along all sight lines. Recall 
that these sight lines emanate from a point which lies in 
the weakest shielding regime; the temperature is ~ 3500K 
where the dissociation rate is calculated, while along the rays 
passing through the dense central core of the halo (for which 
is largest), the temperature drops as low as 300 — 400K. 
The resulting narrowing of the thermal cores in the shield- 
ing gas, along with depopulation of excited rotational states 
from which a large fraction of the dissociating transitions 
occur, drastically reduces the effective column density and 
correspondingly increases fsh- In principle, this effect could 
also contribute to the approximate methods overestimating 
self-shielding as compared to the numerical results; how- 
ever, these changes again are typically only significant for 
the smallest values of the shield factor (sightlines which are 
directed through the dense core), and thus have a small ef- 
fect on the spherically averaged value. 

Finally, we consider the simplified model of a single 
shielding slab with T > Tq. In this case, the larger thermal 
widths in the shielding gas yield greater optical depth in the 
wings of the Voigt profiles of the dissociating transitions, and 
this can outweigh the reduced line centre opacity. This effect 
depends on along the sightline in question, and there 
is a column density-dependent Tmax > To, at which /ah is a 
minimum, above which further increases in the temperature 
lead to weaker shielding. This is not important in the weak- 
est shielding regime, as To here is already very high (see Fig- 
ure [2}, so that few sight lines see a significant increase in the 
temperature. In the strongest shielding regime, the broader 
thermal cores are unimportant because the absorption lines 
already have large damping wings, so that extremely high 
temperatures would be required for this effect to operate. It 
is, however, important in the intermediate regime; as illus- 
trated in the middle panel of Figure [5] This panel shows 
that shielding is always weaker when the temperature is 
held constant along sight lines for /ah, 3D (static) < lO"'^. 
(Here, gas velocities have again been set to zero so as to 
isolate the effect of non-uniform temperature, which cou- 
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pies to that of velocity gradients). However, this cannot be 
entirely attributed to the line-width effect, as temperature 
changes in this regime often also induce stronger shielding 
owing to increased populations of important lines, as men- 
tioned above. Interestingly, the cumulative effect counteracts 
that of Doppler shifts, which are also more important in this 
regime, with the result that "turning on" both temperature 
and velocity gradients causes only a small scatter around 
the original results, as illustrated in the right-hand panel of 
Figure (5] 

3.4 Approximate self— shielding column densities 

The only remaining source of discrepancy between the nu- 
merical and approximate results for fsh are the column 
density approximations, Ai'approx, obtained either from local 
quantities, or from the six-ray method. In order to isolate 
the accuracy of each method described in § 12.3.11 we have 
calculated the dissociation rate numerically, assuming a sin- 
gle slab of shielding gas with A'^h2 = -^approx- The results 
are compared to the exact /ah, 3D, in Figure[71 and the mean 
errors from each method are quantified in Table 2. In the 
three-dimensional calculation, the gas temperature and ve- 
locity are held constant along all sightlines, so as to isolate 
the accuracy of each method for specifying A^approx El- 
Over the full range of radii we consider, the Sobolev 
length proves the most accurate purely local estimate for the 
effective column density. This is perhaps surprising, given 
the results in § 13.21 where we have shown that velocity gra- 
dients do not alter the observed LW flux within these haloes 
significantly. In fact, because of the unimportance of velocity 
gradients for our results, we might expect that the Sobolev 
length should not provide a meaningful length scale, Lchar, 
for self-shielding. It is worth noting that this method is also 
the only one considered here, to the best of our knowledge, 
which has not been implemented in simulations in the con- 
text of self-shielding (though it is commonly used in the 
opposite case, to compute the escape fraction of photons 
traveling outward). Indeed, while it is not obvious how to 
justify the success of this method, its potential usefulness 
in simulations should be highlighted; though the scatter is 
larger than from the six-ray method, the computational ex- 
pense is dramatically reduced. 

In fact. Figure [7] shows that, for /sh < 10~^, all three 
local methods provide surprisingly accurate estimates of the 
column density, given the relatively crude assumptions made 
in approximating the characteristic length scale for shield- 
ing. However, while both Sobolev and Sobolev-like methods 
are also reasonably accurate in the less shielded regime, in 
agreement with SBHIO, we that the Jeans length typically 
underestimates the dissociation rate by an order of magni- 
tude or more at /sh,3D > 10~^, i.e. in regions where the 
number density is low, ntot < 10^ cm~^, and temperature is 
high, T > lO^K as shown by Figure [2] (see also Figure 9 in 
SBH10).~ 

^ We note that it is somewhat awkward to compare the results 
from the Sobolev length to the numerical results assuming a static 
slab — formally the Sobolev length is infinite in static gas. How- 
ever, Figure Vl\ clearly illustrates that the efficacy of this method 
here cannot be attributed to the presence of velocity gradients. 



Table 2. The accuracy of each method for estimating the 
column density (see § 12.311 is quantified by the mean of 
/sh(^approx)//sh 3D (isothermal, static). Ratios are reported sep- 
arately for regions with particle densities above and below 
10'* cm~^, corresponding roughly to a cut— off around /sh^3D ~ 
10~^. Numbers in parentheses indicate the Icr scatter. Boldface 
values highlight the best-fitting Sobolev and Sobolev-like {L'^^^) 
methods. 



ra < 10-^(cm~'') n>10-^(cm-3) 



J cans 


0.50 


(0.38) 


0.96 


(0.17) 




1.7 


(4.3) 


1.1 


(0.3) 


a. r k 


1.3 


(1.3) 


1.0 


(0.3) 


tN 


0.48 


(0.44) 


0.46 


(0.23) 


T min 
-'"Sob 


6.1 


(25) 


2.1 


(0.7) 


h. t' s 
^Sob 


1.4 


(0.9) 


1.1 


(0.2) 


j'k 


0.66 


(0.75) 


0.61 


(0.14) 


t'N 
^Sob 


0.78 


(0.49) 


0.67 


(0.31) 


T min 
^Sob 


3.0 


(6.5) 


1.6 


(0.3) 


six-ray 


0.93 


(0.09) 


0.97 


(0.05) 



Recall that, in this method, the Sobolev length for each sight- 
line is used to calculate the dissociation rate, and the final rate 
is obtained from a simple average over all directions. 

In this case, the rate is calculated given a single column den- 
sity specified by the spherically averaged "Sobolev-like" length. 

Not surprisingly, the six-ray method proves more accu- 
rate than the "local" methods. In fact, in the absence of tem- 
perature and velocity gradients, six sightlines would suffice 
for convergence of the spherically averaged rate. However, 
larger discrepancies occur if the six-ray results are compared 
to the full three-dimensional calculation (this comparison is 
not shown) rather than the static and isothermal case; this 
is simply because the integrated N-n^ is then, in general, not 
equal to the effective self-shielding column density. In this 
case, the average ratio (la scatter) for /sh(six-ray)//sh,3D = 
1.2 (0.1) and 1.1 (0.5) for points with particle densities above 
and below 10* cm~^, respectively (c.f. Table 2). Recall that 
most often this method is implemented by integrating to find 
the total Nyi2 along each of the directions and plugging this 
into equation or ©; therefore, it does not account for 
variations in the temperature and velocity along the rays. 
(An exception is the study bv lYoshida et aL I [2003, in which 
Nn^ is summed along each ray, excluding gas particles with 
relative velocities significantly larger than the local thermal 
velocity.) 

3.5 Angular resolution for /gh in 3D 

In order to determine the required angular resolution for the 
numerical calculations, we performed a convergence test for 
the spherically averaged /sh,3D using 50 points selected from 
the five simulation outputs, spanning the full range of radii 
we consider, 0.1 < -R < lOpc. The mean error in the shield 
factor, /sh,3D, for these points was found to be three (two) 
per cent for an increase from 16 — >■ 25 sightlines (25 49). 
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Figure 7. Comparison of tlie numerical results for the shield factor in three dimensions, /ah, 3D to those obtained from a single slab 
model with (locally) estimated column densities, A^approx- In the former, velocity and temperature are held constant along all sightlines, 
in order to isolate the accuracy of each -/Vapprox method. Clockwise from the upper left, A^approx is specified by the best-fitting Sobolev 
length and Sobolev- like ( "Sobolev-p" ) length (ig^^ and ^g^^,, respectively; see ^ 12.3.111 . and the Jeans length. The lower left panel shows 
the results of the only non-local approach, the "six-ray method," in which the shield factor is calculated for six sight lines using the 
integrated A'^H2 s-i^d average over these gives (six-ray ) . The overall fractional errors of each /sh(Afapprox) arc shown in Table 2. 



In comparison, the mean error for an increase from six to 16 
sightlines was found to be ~ 9 per cent; note that, as in the 
"six-ray" method described above, the six sightlines were 
ahgned parallel with the Cartesian axes, but in this case, 
the full calculation was performed - i.e. including density, 
temperature, and velocity gradients. AH of our calculations 
for /sh,3D are therefore obtained from averaging over 16 sight 
lines, which sample evenly in the azimuthal angle and in the 
cosine of the polar angle. 



An additional convergence test was performed to deter- 
mine the required angular resolution for the mean Sobolev 
and Sobolev-like lengths (L^Sob, Lgob)- For 25 points se- 
lected, again, over the full range of radii considered, the 
mean error in the spherically averaged Sobolev and Sobolev- 
like lengths are 15 and 17 per cent for an increase from 
16 — >■ 25 sightlines, 10 and 11 per cent (25 — >■ 49), and 1 and 
3 per cent (49 — >■ 100), respectively. As a result, 49 sightlines 
were used in all calculations for the Sobolev and Sobolev-like 
methods. (Note that these specific numbers are determined 
by the tiling method we use, which always results in a per- 
fect square for the number of sightlines.) 



4 IMPLICATIONS AND CAVEATS 

4.1 Cooling and collapse in simulated haloes 

In the preceding section, we have shown that common ap- 
proaches to model self-shielding in simulations can intro- 
duce inaccuracies of over an order of magnitude in the H2 
dissociation rate, depending on the level of the shielding and 
on the method being used. Because the cooling properties 
of metal-free gas below 10* K depend sensitively on the H2 
abundance, errors in the dissociation rate can significantly 
affect the thermal and dynamical histories of gas in haloes 
under UV irradiation. In order to quantify such errors, and 
to specify when inaccurate values of /ah will have the greatest 
impact, we show in Figure[S]the H2-cooling and dynamical 
time-scales, given by: 



Tdyn 



'^cool 



16Gp' 



(3/2)nkBT 
I Ah, I ' 



(11) 



(12) 



as functions of radius in the five halo outputs. As in the 
original sim ulations, we employ the H2-cooling rate. Ah,, 
provided by iGaUi fc Palial l| 19981 ): here, G and ks are the 
gravitational and Boltzmann's constant, respectively. 

For all radii considered thus far, 0.1 < 7? < 10 pc, Fig- 
ure[8]shows that the dynamical time exceeds the H2-cooling 
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Figure 8. Upper panel: spherically averaged H2-cooling (solid 
curves) and dynamical (dotted curves) time-scales are sliown for 
the cold haloes; the Hubble time for each output is shown by 
horizontal (dashed) lines, with the same cosmological parameter 
choices as in SBHIO, (f2ni,nA,/i) = (0.279, 0.721,0.701). Thick 
and thin curves denote the results for outputs 1 and 5, respec- 
tively. (The dynamical time-scales are nearly identical for these 
outputs, making the thick and thin curves indistinguishable.) The 
results for output 2 are nearly identical to those of output 1, and 
have been omitted for clarity. Lower panel: same as above for the 
hot haloes, with thick (thin) curves denoting the time-scales for 
output 3 (4). 



Figure 9. Evolution of gas temperature (solid curves and labels 
on left) and /gjj (dashed curves and labels on right) of collapsing 
gas in a halo with T^ir > 10"' K, are obtained from the spherical 
collapse model. Both are shown for runs with intensities, J21, just 
above (thin curves) and below (thick curves) the critical thresh- 
old, Jcrit,2l. The shield factor is obtained from the modified fit- 
ting formula, equation HIOII . and Jcrit 21 = 1-4 x 10^. By contrast, 
if equation ^ is used, Jcrit, 21 = 4.3 X 10'' and the bifurcation 
around threshold flux occurs at /gh 0.1. The incident spectrum 
is modeled as that of a blackbody with T = lO^K and the tem- 
perature of the gas is initialized at the virial temperature of the 
halo. For more on the specifics of the spherical collapse model, 
the reader is referred to SBHIO and references therein. 



time in each of the cold haloes (as expected). More impor- 
tantly, these time-scales are comparable - i.e. within an or- 
der of magnitude - in both hot and cold haloes at 7? < 10 
pc (with the exception of output 3, where the cooling time 
becomes much longer outside the range 2 < _R < 10 pc). 
Because TcooI < Tdyn is required for H2-cooling to impact 
the dynamics of these clouds, it is precisely in regions where 
these time-scales are similar that inaccurate values of the 
dissociation rate may qualitatively change the history of a 
simulated cloud. This will indeed be the case, as long as 
photodissociation and shielding are both non-negligible. For 
example, if the approximated /sh is too large and the H2- 
cooling rate therefore erroneously low, the collapse of over- 
dense clumps in the halo (ultimately leading to formation 
of protostellar objects) may be erroneously delayed. Con- 
versely, if the adopted model overestimates shielding, H2- 
cooling may allow for collapse sooner (on smaller scales) 
than would be possible with the accurate dissociation rate. 

In order to further address how errors in the approx- 
imate /ah may impact the history of a simulated halo, we 
consider the "critical flux," studied by SBHIO. It is well 
known that a sufficiently strong LW flux can suppress H2- 
cooling entirely, keeping the gas in these haloes close to its 
virial temperature, T ~ 10^ K, providing that it remains 
un-enriched in metals. There is a sharp bifurcation in the 
cooling behaviour around the threshold intensity, Jcrit, 21, be- 
low which the gas is able to reach ~ a few hundred Kelvin 
via H2-coohng (outputs 1, 2, and 5), while haloes irradi- 
ated with a super-critical flux will remain H2-poor (outputs 
3 and 4). From the full suite of simulations, SBHIO found 
a critical intensity Jcrit, 21 ~ 10*"^, in the usual units (here 
and throughout this discussion, the results are quoted for an 



incident blackbody spectrum with Tbb ~ 10^ K, relevant for 
direct H2 photodissociation). These results were refined with 
the use of a one-zone spherical collapse model (for details see 
SBHIO, Qmu kai et al. 2008) , and the threshold intensity was 
then found to be Jcrit. 21 ~ 1.2 x lO**. 

To further investigate the implications of inaccurate 
self-shielding models, we have repeated the one-zone cal- 
culations in SBHIO, altering only how /^h is specified. First, 
we adopt the more accurate expression given by DB96, equa- 
tion ([21, rather than their power-law fit, equation (|SJ, em- 
ployed by SBHIO. As in the original study, the H2 column 
density is still specified by the Jeans length (the only vi- 
able option for zero-dimensional models without addition- 
ally computing gradients). The specific intensity, Jlw, is 
then varied iteratively with a Newton-Raphson scheme and 
the critical threshold is found to be Jcrit, 21 = 4.3 x 10'^. 
This reduction in the intensity required to suppress cool- 
ing owes to the temperature dependence of equation (|9]), 
which more accurately models decreased shielding due to 
line-broadening than the power-law. Next, we alter the ex- 
pression for fsh as described in § [ST] (with a = 1-1), and 
find Jcrit, 21 = 1.4 X 10"^. Again, the more accurate fit at high 
temperatures leads to significantly weaker shielding in the 
early stages of collapse. The evolution of both temperature 
and the shield factor as functions of density are shown in 
Figure [9] 

Importantly, because we have employed the Jeans 
length to estimate Nu2, the dissociation rate is still likely 
underestimated by a factor of > 2 at /sh < 10~^ (see the 
upper left panel of Figure [7}. As a result, the critical in- 
tensity quoted here - reduced by a factor of ~ 3 x 2 = 6 
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- is still an upper limit; the full reduction is thus about 
an order of magnitude. Indeed, SBHIO showed that if one 
entirely ignores self-shielding (i.e. the optically thin limit), 
then Jcrit,2i = 4.4 x lO'^ - apparently when self-shielding 
is treated more accurately, it only modestly increases this 
value. 

This reduction in Jcrit,2i illustrates how the self- 
shielding approximations can dramatically impact the ther- 
mal history of a simulated halo. Furthermore, this may also 
have interesting cosmogonical implications, suggesting that 
a larger fraction of haloes than previously thought will see 
a supercritical flux in the spatially fluctuating UV back- 
ground (UVB). It has been suggested (see SBHIO and ref- 
erences therein) that gas irradiated by a supercritical flux 
may avoid fragmentation on small sca les, provided that it 
does not become metal-enriched (see lOmukai et al.l l2008l . 
for the relevant metallicity threshold), and this may provide 
a possible mechanism by which primordial gas could collapse 
directly to form massive black holes, 10*"'' M©. However, 
with the original hig h Jcrit,2i value, it has also been shown 
ijDiikstra et al.|[200 8ll that only one in ~ 10® haloes - only 
those with an unusually bright and close neighbour - will see 
a sufficiently high flux. The direct collapse scenario requires 
that the gas furthermore remain un-poUuted by metals and 
efficiently shed its angular momentum. As a result, it is likely 
that only a small fraction of even these close halo pairs will 
form supermassive black holes. The reduction of the Jcrit,2i 
value will significantly increase the number of candidates for 
objects that avoid H2-cooling and fragmentation, and makes 
this scenario much more viable. 



/sh.HI 



4.2 Internal and external modification of the 
radiation field by HI 

The results presented thus far may be altered by atomic H 
Lyman series absorption either within the haloes, where this 
provides an additional means of shielding H2, or in the IGM, 
where the spectrum of the UVB is modified by the high HI 
optical depth prior to reionization. 

In the first case, suppression of the H2 dissociation rate 
becomes significant if large HI column densities are present, 
A^Hi > 10^^ cm~^. To quantify the impact this has on our 
results for fcdiss, we have repeated the three-dimensional cal- 
culations with the first nine Lyman lines contributing to 
the optical depth of each slab. Note that, while the Lya 
transition (10.2eV) lies outside of the LW wavelength range 
(> 11.2eV), it is nevertheless included, because at large A'hi 
its damping wings protrude into the LW band, and con- 
tribute to LW shielding. 

In the SBHIO haloes, neutral column densities of 



10^ 



cm are typical at R 



10 pc, increasing 



to 10^^~^* cm~^ at radii smaller than a few pc. In these 
densest regions, the dissociation rate is deceased by up to 
a factor of ~ 2.5 in both cold and hot haloes. We find that 
the numerical results for fsii{Nii2, Nhi,T) are well approxi- 
mated by an analytic expression of the form /sh.Ha ^ /sh.Hi, 
where /sh,H2 is the fit for self-shielding in equation 1101 
(with a = 1.1) and /sh.Hi is t he expression provided by 
IWolcott-Green fc HaiinanI (|201ll 'l for HI shielding: 



(1 + a;Hi)' 



X exp (-7 XHi) . 



(13) 



Adopting the same coefHcients as in the original expres- 
sion, l3 = 1.6, 7 = 0.15 and xm = A''hi/2.85 x 10^^ cm"^, 
we find that /sh.H2 x /sh.Hi is accurate to within a factor of 
two in the relevant column density and temperature ranges: 
10" < iVH2 < 10^^ cm-^ 10^^ < iVm < 10^* cm-^ and 
500 < T < 5 X 10^ K. Note that this degree of accuracy is 
perhaps surprising, given that the fit /sh,Hi was originally 
developed for HI shielding of ground state H2 populations; 
for more details on this method for approximating the (non- 
linear) effect of combined shielding by H 2 and HI, the reader 
is referred to IWolcott-Green fc HaimanI (|201ll ). 

The shape of the UVB spectrum is also modified by 
HI in the IGM, as photons traveling over cosmological dis- 
tances will be absorbed once they redshift into resonance 
with a Lyman line. Subsequently, the original photon will 
be replaced, in a radiative cascade, by several lower-energy 
photonfl This resul ts in a characteristi c sawtooth shape 
of the spectrum fsee iHaiman et aLlll997l and HAROO for a 
detailed discussion). To investigate how this modulation of 
the incident flux impacts our results, we have recalculated 
fsh, assuming a single (isothermal and static) slab of gas 
irradiated on one side, with the input spectrum given by 
Ji, = J21X sawtooth modulation (z = 15), adopted from 
HAROO (from the model shown in their Figure 1). We have 
found that the self-shielding behaviour is insensitive to this 
modification - fsh is decreased by ~ 10 per cent over a 
wide range of temperatures (10^ — 10* K), column densi- 
ties (A^H2 < 10'^* cm~^), and irrespective of the assumed 
H2 level populations. While the SBHIO halo outputs are at 
redshifts somewhat lower than z — 15, the magnitude of the 
sawtooth effect decreases over time in the HAROO model, 
due to the redshift -dependent source formation rate; thus, 
the decrease in fsh quoted above is in fact an upper limit|f| 



4.3 Anisotropic LW flux 

In taking the average dissociation rate over all sight lines 
for the 3D calculations, we have thus far assumed isotropy 
of the incident LW flux. This is likely a good approximation 
when the flux is close to the global mean background, and is 
dominated by a large number of distant sources. However, 
if the radiation field is strong, i.e. the halo has one (or a 
few) bright neighbours, then the flux will be anisotropic, 
with fsh for individual sight lines pointing toward the nearby 
source most important. Indeed, the latter case is particularly 
relevant in the context of the critical flux discussed above, 
as supercritical intensiti es, J21 >, 10'^, are li kely to be well 
above the cosmic mean (|Diikstra et al.|[2008l ). 

Let us consider the ramifications of anisotropic incident 



^ In the case of Lya, the original photon may only be lost through 
two-photon decays. However, this line is outside of the wavelength 
range of interest for H2 anyway. 

It is worth noting that the optically thin rate is decreased 
by a factor of ~ 10 compared to that given by the unmodified 
spectrum. However, we are primarily interested here in the self- 
shielding behaviour, rather than in the changes to the model- 
dependent magnitude of the intergalactic UVB. 
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10'= 1Q17 IQIB 10'= lO^O 
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Figure 10. Comparison of the numerically calculated photodis- 
sociation rate coefficient along individual sight lines emanating 
from ~ 100 points at 1 < i? < 10 pc in output 5 (open tri- 
angles) with those obtained when the gas is assumed to be both 
static and isothermal (closed triangles). 

flux for self-shielding models in the context of gas veloc- 
ity and temperature gradients. As shown in Figure [6l the 
shield factor for a given sightline may be altered dramat- 
ically by variations in the gas velocity or temperature in 
that direction. This point is further illustrated in Figure [TOl 
in which /sh for individual sight lines emanating from 100 
points are shown in comparison with those when both ve- 
locity and temperature gradients are artificially switched off. 
(See § I3.2H3.3I for details of how the observed discrepancies 
arise.) Here, we see that the dissociation rate fluctuates by 
«two orders of magnitude in different directions. The rates 
will also be incorrect by up to several orders of magnitude 
if one makes the standard assumption of static and isother- 
mal shielding gas. (On the other hand, the relative accuracy 
of each A'^approx along individual sight lines would not be 
particularly informative here, largely because simulations in 
which the radiation fleld is both spatially and angular re- 
solved typically do not rely on the type of local estimates 
for described above.) 

4.4 Uncertainty in the rovibrational distribution 

In the most general case, one would have to consider the 
non-equilibrium H2 population levels produced during gas- 
phase H2 formation, as well as the time-dependent cascades 
among levels. This, in general, would produce a population 
level distribution that is different than assumed in DB96. 
Absent such a fully time-dependent calculation of the rovi- 
brational level populations, we have thus far made the sim- 
plifying assumption that rotational states within the ground 
vibrational state are populated according to a Boltzmann 
distribution and that the abundance of H2 in higher vibra- 
tional states are negligible. However, because the level pop- 
ulations play an important role, particularly in the temper- 
ature dependence of the dissociation rate (see 13. 11 -^ [33)) . a 
more detailed examination of this assumption is in order. 

It has already been noted that thermal populations are 
not established in higher vibrational levels until particle den- 



sities Merit ~ 10®~^ cm~'^ are reached (see, e.g. Table 1 in 
LPF99), so there is little uncertainty that populations inv > 
are negligible for the SBHIO haloes, at least at the most 
relevant stages of collapse that we considered here. How- 
ever, higher rotational levels within v = Q also thermalize 
at higher critical densities, and because the collisional cross- 
sections depend strongly on temperature, so also do the val- 
ues of ricrit. As an example, consider the S(2) transition (for 
w = 0), for which the critical density at 500, 1000, and 2000 
K, respectively, is Wcrit = 1.2x 10*, 1.7x 10^, 2.5x 10^ cm"^ 
(LPF99). Given that the most rarefied regions of the haloes 
n < 10^ cm~^ are typically at temperatures of ~ several 
thousand Kelvin (see Figure [l]), most likely the rotational 
populations within v = will indeed tend to a Boltzmann 
distribution. None the less, we have repeated our three- 
dimensional calculations with populations in J = 0, 1 only, 
and with a fixed ortho:para ratio of 3:1. The results are com- 
pared to the original /ah, 3D in the left panel of Figure fTTI Be- 
cause restricting the molecules to occupy only J = 0, 1 is, in 
effect, similar to setting the temperature very low (T < 100 
K) in the LTE model, it is not surprising that the shield 
factor is smaller in this case. The clear dip in the results 
above /sh,3D > 5 x 10"* arises simply because the optical 
depth decreases toward hotter regions of these haloes. At the 
smallest values of fsh, the local temperature is only a ~ few 
hundred Kelvin, and as a result, the ground state results do 
not differ dramatically from the original (LTE) calculations. 
In the right panel of Figure fTTI the ground state results are 
also compared to those from equation confirming that 
the DB96 expression is indeed a better fit to the numerical 
results in a rotationally-cold gas. 

4.5 Fluorescent excitation of H2 

All calculations to this point have implicitly assumed that 
each LW photon is permanently removed from the radia- 
tion fleld upon absorption by H2. However, on average only 
10 — 15 per cent of absorption ("pumping") events result 
in dissociation, while the remainder are followed by tran- 
sitions to another bound vibrational state. In a "resonant 
scattering," a single decay returns the molecule directly to 
the initial rovibrational state (w, J = w", J") and the origi- 
nal LW photon is re-emitted. More frequently, UV-pumping 
is followed by a cascade through multiple levels, resulting 
in both infrared fluorescence and emission of a LW pho- 
ton (with different-energy) in the electronic transition. It is 
then possible for this (or the resonantly scattered) photon 
to be reabsorbed, a process which has not been accounted 
for thus far. If the optical depth to the re-emitted photon 
is non-negligible this requires a non-trivial modiflcation in 
our calculations of fcdiss- However, recall that only H2 tran- 
sitions originating in the v = are included; as a result, to 
first order, we need only consider photons emitted in decays 
directly back to the ground vibrational state. The fraction 
of downward t ransit i ons to each bound vibrational level are 
quantified by IShulj I 1978h : they find that ~ 15 per cent 
of all absorption events originating in w = result in de- 
cays directly back to v" = (see their Table 1). Therefore, 
accounting for the optical depth to re-emitted photons rep- 
resents a small correction to our results. 

There is an additional complication, however, if the gas 
is irradiated by a very strong UV fiux, in which the cas- 
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Figure 11. Left panel: comparison of full tliree-dimensional calculations with LTE populations, /sh,3Di to those in which only the ortho 
and para ground states of the molecules are populated, with a fixed ortho:para ratio of 3:1. Right panel: results of the ground-state 
calculation are compared to those from the DB96 fitting formula, equation l[9j. In the latter case, the mean (/gji 37) over all 16 sight lines 
is shown, with velocity and temperature gradients artificially switched off, in order to isolate the accuracy of the fitting formula. 



cade to v" — 0, J" may be interrupted by absorption of LW 
photons from v" > 0. However, only in very strong radi- 
ation fields, J21 > 10^~*, are molecules more likely to be 
"re-pumped" in this way, rather than undergoing radiative 
transitions to lower rovibrational states (|Shulllll978l ). There- 
fore, in the present context of an intergalactic UVB, neglect- 
ing this multiple pumping mechanism is also justified. 



5 SUMMARY AND CONCLUSIONS 

We have shown that the results of existing approximations 
for self-shielding in three-dimensional simulations often 
introduce large inaccuracies in the optically thick H2 
photodissociation rate. In particular, the approximate 
results typically underestimate the numerically calculated 
rate by more than an order of magnitude in low density 
regions, n < 10* cm~^, or where the true shield factor is 
10"'' < /sh,3D < 1. At higher densities, we find that the 
approximate methods are much more accurate, with typical 
errors of ~ 25 per cent. There are a number of factors 
contributing to the discrepancies between the approximate 
and numerical results, which are summarized below: 

(i) The largest source of error in the approximate 
methods is the analytic fit for fsh provided by DB96. 
While this oft-used expression is reasonably accurate at 
low temperatures, we find that it overestimates shielding 
by a large factor at temperatures above a few hundred 
Kelvin (for a gas in LTE). The resulting deviations from 
our numerical results are most apparent at n < 10'* cm"'', 
as these low-density regions typically have not cooled below 
^ 500 - 10^ K in the SBHIO simulations. However, we 
have found that a very simple modification to the DB96 
expression (see equation I10[) improves the agreement with 
our numerical results to within ~ 15 per cent. 

(ii) Nearly all existing approaches to approximate /sh 
are based on a static slab model for the shielding gas, which 
neglects the diminished optical depth due to frequency 
shifts of the H2 resonances in the presence of velocity 
gradients. We find that these frequency shifts do not greatly 



alter the dissociation rate in the SBHIO haloes, largely 
because the gas motions in these simulations are dominated 
by a coherent flow toward the dense core, rather than 
(supersonic) turbulence; therefore, typically only a small 
number of sightlines from a given point see large changes in 
the bulk velocity. The resulting increase in the spherically 
averaged rate is only significant (~ 65 per cent) in the 
range 4 x 10~* < fsh ^ 2 x 10"^, and elsewhere is generally 
negligible. 

(iii) We find that gas temperature gradients can alter 
the self-shielding behavior dramatically via changes in 
both the thermal line widths and in the rovibrational level 
populations. In this case, the dissociation rate may be 
either increased or decreased, depending on the sign of the 
temperature gradient along each sightline. However, the 
effect again is smaller for the spherically averaged rate, 
and tends to make shielding stronger in the haloes we 
analyze, largely counteracting any changes due to frequency 
shifts. Taken together, temperature and velocity gradients 
therefore typically only introduce a factor of < 2 scatter in 
spherically averaged rate. 

(iv) Finally, we have evaluated several approaches to esti- 
mate the shielding column density in simulations (and one- 
zone models). The most common of these is based on the 
assumption that the characteristic length scale for shield- 
ing is of order the local Jeans length. In agreement with 
SBHIO, we find that this method is very accurate at den- 
sities n > 10*"^ cm~^, but underestimates the optically 
thick rate by > an order of magnitude at lower densities. 
However, we show that two less commonly used methods 
provide reasonably accurate estimates of the shielding col- 
umn at all densities we consider, and are computationally 
inexpensive, relying again only on local properties of the gas. 
These are based on the Sobolev length and a variation of the 
Sobolev length based on the density - rather than velocity 
- gradient; both yield more accurate results for f^h than the 
Jeans length method in low-density regions, with essentially 
unbiased scatter around the true value. We also show that 
a "six-ray" method, based on integrating the column den- 
sity in only six directions is extremely accurate, deviating 
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from the "exact results" only due to the effects of tempera- 
ture and velocity gradients. However, this non-local method 
comes with a larger computational expense. 

In addition to the factors enumerated above, HI shield- 
ing of H2 also causes the "true" dissociation rate to deviate 
from the results of approximate treatments, which most of- 
ten neglect this additional shielding. We find that this effect 
can be w ell modeled by a simple analytic prescription pro- 
vided by IWolcott-Green fc HaimanI l|201ll '). The simple fit- 
ting formulae we provide can be trivially incorporated into 
future three-dimensional simulations, to improve the speed 
and accuracy of calculations of the H2-photodissociation 
rate. 

Finally, using the same (one-zone) spherical-collapse 
model employed by SBHIO, we show that the critical LW 
flux Jcrit required to keep Tvir > lO^K haloes H2-poor is 
reduced by about an order of magnitude as a result of im- 
proved accuracy in the shielding factor. This serves to illus- 
trate an important conclusion of our study: the cooling prop- 
erties - and thereby the dynamical history - of (metal-free) 
gas in simulated haloes depend sensitively on the adopted 
self-shielding model. In particular, the reduction in Jcrit im- 
plies that H2-cooling is suppressed in many more of these 
haloes, thus increasing the potential sites for direct forma- 
tion of SMBHs in the early universe. 
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